# Install packages
# library(remotes)
# remotes::install_github('ropensci/osmdata')
# Data manipulation
# library(tibble)
library(dplyr)
# library(reshape)
# Data visualisation
library(ggplot2)
# Working with spatial data
# library(sf)
# library(leaflet)
library(raster)
library(rgdal)
# library(proj4)
# library(scales)
# Working with OSM data
# library(osmdata)
# library(osmextract)
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area
HARV_dsmCrop_info <- capture.output(
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
)
We are going to work with a Digital Surface Model (DSM) which is in the GeoTIFF format.
DSM_HARV <-
raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
DSM_HARV
## class : RasterLayer
## dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_dsmCrop.tif
## names : HARV_dsmCrop
## values : 305.07, 416.07 (min, max)
summary(DSM_HARV)
## HARV_dsmCrop
## Min. 305.5500
## 1st Qu. 345.6500
## Median 359.6450
## 3rd Qu. 374.2825
## Max. 413.9000
## NA's 0.0000
summary(DSM_HARV, maxsamp = ncell(DSM_HARV))
## HARV_dsmCrop
## Min. 305.07
## 1st Qu. 345.59
## Median 359.67
## 3rd Qu. 374.28
## Max. 416.07
## NA's 0.00
DSM_HARV_df <- as.data.frame(DSM_HARV, xy = TRUE)
str(DSM_HARV_df)
## 'data.frame': 2319799 obs. of 3 variables:
## $ x : num 731454 731454 731456 731456 731458 ...
## $ y : num 4713838 4713838 4713838 4713838 4713838 ...
## $ HARV_dsmCrop: num 409 408 407 407 409 ...
ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_viridis_c() +
coord_quickmap()
plot(DSM_HARV)
But what units are these?
# CRS have been explained on Day 1
crs(DSM_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
## BBOX[0,-78,84,-72]],
## ID["EPSG",32618]]
minValue(DSM_HARV)
## [1] 305.07
maxValue(DSM_HARV)
## [1] 416.07
If Min and Max values haven’t been calculated, you can set them with
the raster::setMinMax() function.
# Maybe skip this
DSM_HARV <- raster::setMinMax(DSM_HARV)
To see how many bands a raster dataset has, use the
raster::nlayers() function.
nlayers(DSM_HARV)
## [1] 1
This dataset has only 1 band. We will discuss multi-band raster data in a later episode.
In raster data, pixels with no value are represented with
NoDataValue. That is usually the case with edge values
where data is missing from the rectangular region of the raster. Use
raster::GDALinfo() to check the metadata for missing data
values.
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area
Bad values are usually miscalculations that result in values out of a predefined range.
# This code chunk is hidden in the lesson material
A histogram can be used to inspect the distribution of raster values
visually. It can show if there are values above the max or below the min
of the expected range. We can create a histogram with the ggplot2
function geom_histogram().
ggplot() +
geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop))
Adjust the level of desired detail by setting the number of bins.
ggplot() +
geom_histogram(data = DSM_HARV_df, aes(HARV_dsmCrop), bins = 40)
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 -0.7136298 0.9999997 0.3125525 0.4812939
## Metadata:
## AREA_OR_POINT=Area
DSM_HARV_df <- DSM_HARV_df %>%
mutate(fct_elevation = cut(HARV_dsmCrop, breaks = 3))
ggplot() +
geom_bar(data = DSM_HARV_df, aes(fct_elevation))
To see the cutoff values:
unique(DSM_HARV_df$fct_elevation)
## [1] (379,416] (342,379] (305,342]
## Levels: (305,342] (342,379] (379,416]
To show count number of pixels in each group:
DSM_HARV_df %>%
group_by(fct_elevation) %>%
count()
## # A tibble: 3 × 2
## # Groups: fct_elevation [3]
## fct_elevation n
## <fct> <int>
## 1 (305,342] 418891
## 2 (342,379] 1530073
## 3 (379,416] 370835
To customize cutoff values:
custom_bins <- c(300, 350, 400, 450)
DSM_HARV_df <- DSM_HARV_df %>%
mutate(fct_elevation_2 = cut(HARV_dsmCrop, breaks = custom_bins))
unique(DSM_HARV_df$fct_elevation_2)
## [1] (400,450] (350,400] (300,350]
## Levels: (300,350] (350,400] (400,450]
ggplot() +
geom_bar(data = DSM_HARV_df, aes(fct_elevation_2))
DSM_HARV_df %>%
group_by(fct_elevation_2) %>%
count()
## # A tibble: 3 × 2
## # Groups: fct_elevation_2 [3]
## fct_elevation_2 n
## <fct> <int>
## 1 (300,350] 741815
## 2 (350,400] 1567316
## 3 (400,450] 10668
ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = fct_elevation_2)) +
coord_quickmap()
Customising colors:
terrain.colors(3)
## [1] "#00A600" "#ECB176" "#F2F2F2"
ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y,
fill = fct_elevation_2)) +
scale_fill_manual(values = terrain.colors(3)) +
coord_quickmap()
my_col <- terrain.colors(3)
ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y,
fill = fct_elevation_2)) +
scale_fill_manual(values = my_col, name = "Elevation") +
coord_quickmap()
DSM_HARV_df <- DSM_HARV_df %>%
mutate(fct_elevation_6 = cut(HARV_dsmCrop, breaks = 6))
unique(DSM_HARV_df$fct_elevation_6)
## [1] (398,416] (379,398] (361,379] (342,361] (324,342] (305,324]
## Levels: (305,324] (324,342] (342,361] (361,379] (379,398] (398,416]
my_col <- terrain.colors(6)
ggplot() +
geom_raster(data = DSM_HARV_df, aes(x = x, y = y,
fill = fct_elevation_6)) +
scale_fill_manual(values = my_col, name = "Elevation") +
coord_quickmap() +
xlab("X") +
ylab("Y") +
labs(title = "Elevation Classes of the Harvard Forest Digital Surface Model (DSM)")
DSM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
DSM_hill_HARV
## class : RasterLayer
## dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_DSMhill.tif
## names : HARV_DSMhill
## values : -0.7136298, 0.9999997 (min, max)
DSM_hill_HARV_df <- as.data.frame(DSM_hill_HARV, xy = TRUE)
str(DSM_hill_HARV_df)
## 'data.frame': 2319799 obs. of 3 variables:
## $ x : num 731454 731454 731456 731456 731458 ...
## $ y : num 4713838 4713838 4713838 4713838 4713838 ...
## $ HARV_DSMhill: num NA NA NA NA NA NA NA NA NA NA ...
ggplot() +
geom_raster(data = DSM_hill_HARV_df,
aes(x = x, y = y, alpha = HARV_DSMhill)) +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
coord_quickmap()
ggplot() +
geom_raster(data = DSM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dsmCrop)) +
geom_raster(data = DSM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DSMhill)) +
scale_fill_viridis_c() +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
ggtitle("Elevation with hillshade") +
coord_quickmap()
# CREATE DSM MAPS
# import DSM data
DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
# convert to a df for plotting
DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)
# import DSM hillshade
DSM_hill_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmHill.tif")
# convert to a df for plotting
DSM_hill_SJER_df <- as.data.frame(DSM_hill_SJER, xy = TRUE)
# Build Plot
ggplot() +
geom_raster(data = DSM_SJER_df ,
aes(x = x, y = y,
fill = SJER_dsmCrop,
alpha = 0.8)
) +
geom_raster(data = DSM_hill_SJER_df,
aes(x = x, y = y,
alpha = SJER_dsmHill)
) +
scale_fill_viridis_c() +
guides(fill = guide_colorbar()) +
scale_alpha(range = c(0.4, 0.7), guide = "none") +
# remove grey background and grid lines
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
xlab("UTM Easting Coordinate (m)") +
ylab("UTM Northing Coordinate (m)") +
ggtitle("DSM with Hillshade") +
coord_quickmap()
# CREATE DTM MAP
# import DTM
DTM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmCrop.tif")
DTM_SJER_df <- as.data.frame(DTM_SJER, xy = TRUE)
# DTM Hillshade
DTM_hill_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_dtmHill.tif")
DTM_hill_SJER_df <- as.data.frame(DTM_hill_SJER, xy = TRUE)
ggplot() +
geom_raster(data = DTM_SJER_df ,
aes(x = x, y = y,
fill = SJER_dtmCrop,
alpha = 2.0)
) +
geom_raster(data = DTM_hill_SJER_df,
aes(x = x, y = y,
alpha = SJER_dtmHill)
) +
scale_fill_viridis_c() +
guides(fill = guide_colorbar()) +
scale_alpha(range = c(0.4, 0.7), guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
ggtitle("DTM with Hillshade") +
coord_quickmap()
What happens when maps don’t line up? That is usually a sign that layers are in different coordinate reference systems (CRS).
DTM_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
DTM_hill_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")
DTM_HARV_df <- as.data.frame(DTM_HARV, xy = TRUE)
DTM_hill_HARV_df <- as.data.frame(DTM_hill_HARV, xy = TRUE)
ggplot() +
geom_raster(data = DTM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dtmCrop)) +
geom_raster(data = DTM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DTMhill_WGS84)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
ggplot() +
geom_raster(data = DTM_HARV_df,
aes(x = x, y = y,
fill = HARV_dtmCrop)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
ggplot() +
geom_raster(data = DTM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DTMhill_WGS84)) +
coord_quickmap()
### Challenge
crs(DTM_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
## BBOX[0,-78,84,-72]],
## ID["EPSG",32618]]
crs(DTM_hill_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## REMARK["Axis order reversed compared to EPSG:4326"]]
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = crs(DTM_HARV))
crs(DTM_hill_UTMZ18N_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
## BBOX[0,-78,84,-72]],
## ID["EPSG",32618]]
crs(DTM_hill_HARV)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## REMARK["Axis order reversed compared to EPSG:4326"]]
extent(DTM_hill_UTMZ18N_HARV)
## class : Extent
## xmin : 731397.3
## xmax : 733205.3
## ymin : 4712403
## ymax : 4713907
extent(DTM_hill_HARV)
## class : Extent
## xmin : -72.18192
## xmax : -72.16061
## ymin : 42.52941
## ymax : 42.54234
res(DTM_hill_UTMZ18N_HARV)
## [1] 1.000 0.998
res(DTM_HARV)
## [1] 1 1
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = crs(DTM_HARV),
res = res(DTM_HARV))
res(DTM_hill_UTMZ18N_HARV)
## [1] 1 1
res(DTM_HARV)
## [1] 1 1
DTM_hill_HARV_2_df <- as.data.frame(DTM_hill_UTMZ18N_HARV, xy = TRUE)
ggplot() +
geom_raster(data = DTM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dtmCrop)) +
geom_raster(data = DTM_hill_HARV_2_df,
aes(x = x, y = y,
alpha = HARV_DTMhill_WGS84)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
# import DSM
DSM_SJER <- raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
# import DSM hillshade
DSM_hill_SJER_WGS <-
raster("data/NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")
# reproject raster
DTM_hill_UTMZ18N_SJER <- projectRaster(DSM_hill_SJER_WGS,
crs = crs(DSM_SJER),
res = 1)
# convert to data.frames
DSM_SJER_df <- as.data.frame(DSM_SJER, xy = TRUE)
DSM_hill_SJER_df <- as.data.frame(DTM_hill_UTMZ18N_SJER, xy = TRUE)
ggplot() +
geom_raster(data = DSM_hill_SJER_df,
aes(x = x, y = y,
alpha = SJER_DSMhill_WGS84)
) +
geom_raster(data = DSM_SJER_df,
aes(x = x, y = y,
fill = SJER_dsmCrop,
alpha=0.8)
) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 304.56 389.82 344.8979 15.86147
## Metadata:
## AREA_OR_POINT=Area
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
## rows 1367
## columns 1697
## bands 1
## lower left origin.x 731453
## lower left origin.y 4712471
## res.x 1
## res.y 1
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 1697
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 305.07 416.07 359.8531 17.83169
## Metadata:
## AREA_OR_POINT=Area
ggplot() +
geom_raster(data = DTM_HARV_df ,
aes(x = x, y = y, fill = HARV_dtmCrop)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
ggplot() +
geom_raster(data = DSM_HARV_df ,
aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
CHM_HARV <- DSM_HARV - DTM_HARV
CHM_HARV_df <- as.data.frame(CHM_HARV, xy = TRUE)
ggplot() +
geom_raster(data = CHM_HARV_df ,
aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_quickmap()
ggplot(CHM_HARV_df) +
geom_histogram(aes(layer))
custom_bins <- c(0, 10, 20, 30, 40)
CHM_HARV_df <- CHM_HARV_df %>%
mutate(canopy_discrete = cut(layer, breaks = custom_bins))
ggplot() +
geom_raster(data = CHM_HARV_df , aes(x = x, y = y,
fill = canopy_discrete)) +
scale_fill_manual(values = terrain.colors(4)) +
coord_quickmap()
CHM_ov_HARV <- overlay(DSM_HARV,
DTM_HARV,
fun = function(r1, r2) { return( r1 - r2) })
CHM_ov_HARV_df <- as.data.frame(CHM_ov_HARV, xy = TRUE)
ggplot() +
geom_raster(data = CHM_ov_HARV_df,
aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_quickmap()
writeRaster(CHM_ov_HARV, "CHM_HARV.tiff",
format="GTiff",
overwrite=TRUE,
NAflag=-9999)
CHM_ov_SJER <- overlay(DSM_SJER,
DTM_SJER,
fun = function(r1, r2){ return(r1 - r2) })
CHM_ov_SJER_df <- as.data.frame(CHM_ov_SJER, xy = TRUE)
ggplot(CHM_ov_SJER_df) +
geom_histogram(aes(layer))
ggplot() +
geom_raster(data = CHM_ov_SJER_df,
aes(x = x, y = y,
fill = layer)) +
scale_fill_gradientn(name = "Canopy Height",
colors = terrain.colors(10)) +
coord_quickmap()
writeRaster(CHM_ov_SJER, "chm_ov_SJER.tiff",
format = "GTiff",
overwrite = TRUE,
NAflag = -9999)
ggplot(CHM_HARV_df) +
geom_histogram(aes(layer))
ggplot(CHM_ov_SJER_df) +
geom_histogram(aes(layer))
##Getting Started with Multi-Band Data in R
The raster() function only reads in the first band, in
this case the red band of an RGB raster
RGB_band1_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
RGB_band1_HARV_df <- as.data.frame(RGB_band1_HARV, xy = TRUE)
ggplot() +
geom_raster(data = RGB_band1_HARV_df,
aes(x = x, y = y, alpha = HARV_RGB_Ortho)) +
coord_quickmap()
RGB_band1_HARV
## class : RasterLayer
## band : 1 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho
## values : 0, 255 (min, max)
nbands(RGB_band1_HARV)
## [1] 3
To import the green band:
RGB_band2_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif", band = 2)
RGB_band2_HARV_df <- as.data.frame(RGB_band2_HARV, xy = TRUE)
ggplot() +
geom_raster(data = RGB_band2_HARV_df,
aes(x = x, y = y, alpha = HARV_RGB_Ortho)) +
coord_equal()
The stack() function brings in all bands
RGB_stack_HARV <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
RGB_stack_HARV
## class : RasterStack
## dimensions : 2317, 3073, 7120141, 3 (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : HARV_RGB_Ortho.1, HARV_RGB_Ortho.2, HARV_RGB_Ortho.3
## min values : 0, 0, 0
## max values : 255, 255, 255
RGB_stack_HARV@layers
## [[1]]
## class : RasterLayer
## band : 1 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.1
## values : 0, 255 (min, max)
##
##
## [[2]]
## class : RasterLayer
## band : 2 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.2
## values : 0, 255 (min, max)
##
##
## [[3]]
## class : RasterLayer
## band : 3 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.3
## values : 0, 255 (min, max)
RGB_stack_HARV[[2]]
## class : RasterLayer
## band : 2 (of 3 bands)
## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : HARV_RGB_Ortho.tif
## names : HARV_RGB_Ortho.2
## values : 0, 255 (min, max)
RGB_stack_HARV_df <- as.data.frame(RGB_stack_HARV, xy = TRUE)
str(RGB_stack_HARV_df)
## 'data.frame': 7120141 obs. of 5 variables:
## $ x : num 731999 731999 731999 731999 732000 ...
## $ y : num 4713535 4713535 4713535 4713535 4713535 ...
## $ HARV_RGB_Ortho.1: num 0 2 6 0 16 0 0 6 1 5 ...
## $ HARV_RGB_Ortho.2: num 1 0 9 0 5 0 4 2 1 0 ...
## $ HARV_RGB_Ortho.3: num 0 10 1 0 17 0 4 0 0 7 ...
ggplot() +
geom_histogram(data = RGB_stack_HARV_df, aes(HARV_RGB_Ortho.1))
ggplot() +
geom_raster(data = RGB_stack_HARV_df,
aes(x = x, y = y, alpha = HARV_RGB_Ortho.2)) +
coord_quickmap()
plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3)
Which is the same with terra
terra::plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3)
plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3,
scale = 800,
stretch = "lin")
plotRGB(RGB_stack_HARV,
r = 1, g = 2, b = 3,
scale = 800,
stretch = "hist")
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")
## rows 2317
## columns 3073
## bands 3
## lower left origin.x 731998.5
## lower left origin.y 4712956
## res.x 0.25
## res.y 0.25
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -9999 1 3073
## 2 Float64 TRUE -9999 1 3073
## 3 Float64 TRUE -9999 1 3073
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 0 255 107.83651 30.01918
## 2 0 255 130.09595 32.00168
## 3 0 255 95.75979 16.57704
## Metadata:
## AREA_OR_POINT=Area
HARV_NA <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")
plotRGB(HARV_NA,
r = 1, g = 2, b = 3)
GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
## rows 2317
## columns 3073
## bands 3
## lower left origin.x 731998.5
## lower left origin.y 4712956
## res.x 0.25
## res.y 0.25
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## file data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -1.7e+308 1 3073
## 2 Float64 TRUE -1.7e+308 1 3073
## 3 Float64 TRUE -1.7e+308 1 3073
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 0 255 NaN NaN
## 2 0 255 NaN NaN
## 3 0 255 NaN NaN
## Metadata:
## AREA_OR_POINT=Area
object.size(RGB_stack_HARV)
## 51296 bytes
RGB_brick_HARV <- brick(RGB_stack_HARV)
object.size(RGB_brick_HARV)
## 170898912 bytes
plotRGB(RGB_brick_HARV)
methods(class=class(RGB_stack_HARV))
## [1] ! != [
## [4] [[ [[<- [<-
## [7] %in% == $
## [10] $<- addLayer adjacent
## [13] aggregate all.equal animate
## [16] approxNA area Arith
## [19] as.array as.character as.data.frame
## [22] as.integer as.list as.logical
## [25] as.matrix as.vector atan2
## [28] bbox blockSize boxplot
## [31] brick calc cellFromRowCol
## [34] cellFromRowColCombine cellFromXY cellStats
## [37] clamp click coerce
## [40] colFromCell colFromX colSums
## [43] Compare coordinates corLocal
## [46] couldBeLonLat cover crop
## [49] crosstab crs<- cut
## [52] cv density dim
## [55] dim<- disaggregate dropLayer
## [58] extend extent extract
## [61] flip freq getValues
## [64] getValuesBlock getValuesFocal hasValues
## [67] head hist image
## [70] init inMemory interpolate
## [73] intersect is.factor is.finite
## [76] is.infinite is.na is.nan
## [79] isLonLat KML labels
## [82] length levels levels<-
## [85] log Logic mask
## [88] match Math Math2
## [91] maxValue mean merge
## [94] metadata minValue modal
## [97] mosaic names names<-
## [100] ncell ncol ncol<-
## [103] nlayers nrow nrow<-
## [106] origin origin<- overlay
## [109] pairs persp plot
## [112] plotRGB predict print
## [115] proj4string proj4string<- quantile
## [118] raster rasterize ratify
## [121] readAll readStart readStop
## [124] reclassify rectify res
## [127] res<- resample rotate
## [130] rowColFromCell rowFromCell rowFromY
## [133] rowSums sampleRandom sampleRegular
## [136] scale select setMinMax
## [139] setValues shift show
## [142] spplot stack stackSelect
## [145] stretch subs subset
## [148] Summary summary t
## [151] tail text trim
## [154] unique unstack values
## [157] values<- weighted.mean which.max
## [160] which.min whiches.max whiches.min
## [163] wkt writeRaster xFromCell
## [166] xFromCol xmax xmax<-
## [169] xmin xmin<- xres
## [172] xyFromCell yFromCell yFromRow
## [175] ymax ymax<- ymin
## [178] ymin<- yres zonal
## [181] zoom
## see '?methods' for accessing help and source code
methods(class=class(RGB_stack_HARV[1]))
## [1] [ [<- anyDuplicated as_tibble as.data.frame
## [6] as.raster bind boxplot brick cellFromXY
## [11] coerce coordinates determinant distance duplicated
## [16] edit extent extract head initialize
## [21] isSymmetric Math Math2 Ops raster
## [26] rasterize relist subset summary surfaceArea
## [31] tail trim unique values<- weighted.mean
## [36] writeValues
## see '?methods' for accessing help and source code